setwd("/home1/30/jc227089/Quoll_model/Random_outputs")

# vector of observations: popsize and sexratio at Pob and Astell islands
obs<-c(766, 818, 470, 300, 64/(18), 87/(55), 46/(2), 30/(4), 3228, 4820, 2590, 2890, 134/(8), 167/(24), 84/1, 103/5)

# function to standardize a vector (mean=0, sd=1)
stdz<-function(vec){
	(vec-mean(vec))/sd(vec)	
}

gold<-c() # place to put the best models
for (i in 1:100){
	fname=paste("ABCSample", i, ".RData", sep="")
	load(fname)
	ID<-1:length(out[,1])
	out<-cbind(out, ID)
	temp<-subset(out, is.na(out[,"pob.sr7"])==F & is.na(out[,"as.sr7"])==F) #find extinctions and remove
	temp<-subset(temp, apply(temp[,c(12:15, 20:23)], 1, max)!=Inf & apply(temp[,c(12:15, 20:23)], 1, min)!=0) #find sex ratios = inf or 0 and remove
	temp2<-temp[,8:23] #gather model data
	qtile<-0.005*5000/length(temp[,1]) #work out quantile to sample
	if (is.na(qtile)==T) {
		rm(out)
		next	
	}
	temp2<-sweep(temp2, 2, obs) # calculate sum of squares
	temp2<-temp2^2 
	ss.pop.pb<-apply(temp2[,c(1:4)], 1, sum) 
	ss.pop.as<-apply(temp2[,c(9:12)], 1, sum) 
	ss.pop<-stdz(ss.pop.pb)+stdz(ss.pop.as) #sum standardised deviations
	test<-ss.pop
	quant.pop<-quantile(test, qtile) #find quantiles
	keepers<-which(test<=quant.pop) 
	keepers<-match(temp[keepers,"ID"], out[,"ID"])
	gold<-rbind(gold, out[keepers,]) #keep these rows and rbind to gold
	rm(out)	
}

save(gold, file="Kept_sims_poponlyhalfpcED.RData")